In [28]:
#核心论文代码验证
using Plots, ComplexPhasePortrait, ApproxFun, SpecialFunctions
gr();

u is in complex plane.

λ02=(1+uM0)2u2
λ12=C12(1+uM0)2u2
The principal branches of the square roots are used so that the branch cuts of λ0± go from u0± to and the branch cuts of λ1± go from u1± to where we have defined u0±=1/(M0±1) and u1±=C1/(M1C1±1).

In [2]:
#u is in complex plane,for lanta_0,it has branch cuts
M0=0.3
phaseplot(-3..3, -3..3, z -> sqrt((1-z*M0)^2-z^2))
Out[2]:
-3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3

The fourier transform of the hard-wall boundary condtion on the centre body is simply β/r(h,u)=0满足. To satisfy this boundary condition, we choose to write the solution of Bessel eqution as follows:

β(r,u)={A(u))Hm(1)(λ0wr) if r>1B(u)[Ym(λ1wh)Jm(λ1wr)Jm(λ1wh)Ym(λ1wr)] if h<r<1

where Jm and Ym are the Bessel and Neumann functions of order m, and Hm(1)=Jm+iYm is the Hankel function of the first type.

参考:关于bessel的导数:https://www.boost.org/doc/libs/1_57_0/libs/math/doc/html/math_toolkit/bessel/bessel_derivatives.html http://functions.wolfram.com/Bessel-TypeFunctions/BesselJ/20/ShowAll.html

Jm=Jm1Jm+12Ym=Ym1Ym+12

In [3]:
#首先让我们大致看一下bessel函数的曲线形式
#与https://dlmf.nist.gov/10.3一致

𝐽𝑚=(m,x) -> besselj(m, x)
𝑌𝑚=(m,x) -> bessely(m, x)
x = range(0.01,stop=10,length=1000)
plot(x, 𝐽𝑚.(0,x);label="J_0",ylims = (-1,1))
plot!(x, 𝐽𝑚.(1,x);label="J_1")
plot!(x, 𝑌𝑚.(0,x);label="Y_0")
plot!(x, 𝑌𝑚.(1,x);label="Y_1")
Out[3]:
0.0 2.5 5.0 7.5 10.0 -1.0 -0.5 0.0 0.5 1.0 J_0 J_1 Y_0 Y_1
In [4]:
m=1
𝐽′𝑚=(m,x) -> 1/2*(besselj(m-1, x)-besselj(m+1, x))
𝑌′𝑚=(m,x) -> 1/2*(bessely(m-1, x)-bessely(m+1, x))
x = range(0.01,stop=10,length=1000)
plot(x, 𝐽′𝑚.(m,x);label="J'_1")
plot!(x, 𝑌′𝑚.(m,x);label="Y'_1",ylims = (-1,1))
Out[4]:
0.0 2.5 5.0 7.5 10.0 -1.0 -0.5 0.0 0.5 1.0 J'_1 Y'_1
In [5]:
f = z -> hankelh1(3,z)
# set up plotting grid
z = range(0.01 ,stop=100,  length=1000)
plot(z, real(f.(z)))
Out[5]:
0 25 50 75 100 -0.2 0.0 0.2 0.4 y1
In [74]:
#验证汉克函数
𝐻𝑚=(m,z) -> besselj(m, z)+im*bessely(m, z)
@show 𝐻𝑚(2, 1+im)
@show hankelh1(2,1+im)
phaseplot(-3..3, -3..3, z -> hankelh1(2,z))
#验证汉克函数在r趋近与无穷大时,叙述为正和负的区别
@show hankelh1(1,-100im)
@show hankelh1(1,+100im)
𝐻𝑚(2, 1 + im) = -0.5357570706365329 - 0.22597037902118686im
hankelh1(2, 1 + im) = -0.5357570706365328 - 0.22597037902118686im
hankelh1(1, -100im) = 1.30837515418464e26 - 2.1367387806763253e42im
hankelh1(1, 100im) = -2.9792874198947446e-45 - 1.8242874012570376e-61im
Out[74]:
-2.9792874198947446e-45 - 1.8242874012570376e-61im

必须要指出的是,hankel函数用Jm+imYm近似的方法只能针对与实数轴,对于过大的正虚数,imag》20im,则将不成立! 原因是,hankelh1(1, 20im)-->0,而besselj(1, 20im)-->∞&bessely(1, 20im)-->∞,因此造成很大的数值误差. 举例说明:

In [259]:
#验证hankel函数的导数
𝐻1𝑚=(m,z) -> hankelh1(m, z)
H1𝑚=(m,z) -> besselj(m, z)+im*bessely(m, z)
𝐻′1𝑚1=(m,z) ->hankelh1(m-1, z)-m/z*hankelh1(m, z)
H′1𝑚1=(m,z) ->besselj(m-1, z)+im*bessely(m-1, z)-m/z*(besselj(m, z)+im*bessely(m, z))
#H′1𝑚2=(m,z) ->-(besselj(m+1, z)+im*bessely(m+1, z))+m/z*(besselj(m, z)+im*bessely(m, z))
#𝐻′1𝑚2=(m,z) ->-hankelh1(m+1, z)+m/z*hankelh1(m, z)

@show besselj(1,2+20*im)
@show bessely(1,2+20*im)
@show 𝐻1𝑚(1,2+20*im)
@show H1𝑚(1,2+20*im)
@show besselj(1,2-20*im)
@show bessely(1,2-20*im)
@show 𝐻1𝑚(1,2-20*im)
@show H1𝑚(1,2-20*im)

zr=range(-40,stop=40,length=1000)
zi=ones(1000,1)*im*20
z=zr+zi
r_j = z -> real(besselj(1,z))
i_j = z -> imag(besselj(1,z))
r_y = z -> real(bessely(1,z))
i_y = z -> imag(bessely(1,z))
r_h=z -> real(hankelh1(1,z))
i_h=z -> imag(hankelh1(1,z))
#plot(zr,r_j.(z);label="r-j")
#plot!(zr,i_j.(z);label="i-j")
#plot!(zr,r_y.(z);label="r-y")
#plot!(zr,i_y.(z);label="i-y")
plot(zr,r_h.(z);label="T-r-h")
plot!(zr,i_h.(z);label="T-i-h")

#plot!(zr,r_j.(z)+i_y.(z);label="F-r-h")
#plot!(zr,i_j.(z)+r_y.(z);label="F-i-h")
besselj(1, 2 + 20im) = 3.931553698919195e7 - 1.576352471413e7im
bessely(1, 2 + 20im) = 1.5763524714130001e7 + 3.931553698919195e7im
𝐻1𝑚(1, 2 + 20im) = 1.727569618545747e-10 - 3.311815065025103e-10im
H1𝑚(1, 2 + 20im) = 0.0 + 1.862645149230957e-9im
besselj(1, 2 - 20im) = 3.931553698919195e7 + 1.576352471413e7im
bessely(1, 2 - 20im) = 1.5763524714130001e7 - 3.931553698919195e7im
𝐻1𝑚(1, 2 - 20im) = 7.86310739783839e7 + 3.1527049428260002e7im
H1𝑚(1, 2 - 20im) = 7.86310739783839e7 + 3.152704942826e7im
Out[259]:
-40 -20 0 20 40 - 3×10 - 10 - 2×10 - 10 - 1×10 - 10 0 1×10 - 10 2×10 - 10 3×10 - 10 T-r-h T-i-h
In [6]:
M0=0.3
m=2
w=1+im
λ0=(u)->sqrt((1-u*M0)^2-u^2)
𝐻1𝑚=(u) -> besselj(m,u)+im*bessely(m,u)
𝐻2𝑚=(u) -> hankelh1(m,u)
𝐻1𝑚_λ0w=(u) -> besselj(m, λ0(u)*w)+im*bessely(m, λ0(u)*w)
𝐻2𝑚_λ0w=(u) -> hankelh1(m, λ0(u)*w)
𝐻′1𝑚1_λ0w=(u) ->hankelh1(m-1, λ0(u)*w)-m/(λ0(u)*w)*hankelh1(m, λ0(u)*w)
H′1𝑚1_λ0w=(u) ->besselj(m-1, λ0(u)*w)+im*bessely(m-1, λ0(u)*w)-m/(λ0(u)*w)*(besselj(m, λ0(u)*w)+im*bessely(m, λ0(u)*w))

phaseplot(-20..20, -40..40, 𝐻′1𝑚1_λ0w)
phaseplot(-20..20, -40..40, H′1𝑚1_λ0w)
Out[6]:
-20 -10 0 10 20 -40 -20 0 20 40
In [255]:
m=2
w=1+im
λ0=(u)->sqrt((1-u*M0)^2-u^2)
𝐻1𝑚_λ0w=(u) -> besselj(m, λ0(u)*w)+im*bessely(m, λ0(u)*w)
𝐻2𝑚_λ0w=(u) -> hankelh1(m, λ0(u)*w)
phaseplot(-20..20, -20..20, 𝐻2𝑚_λ0w)
phaseplot!(-20..20, -20..20, 𝐻1𝑚_λ0w)
Out[255]:
-20 -10 0 10 20 -20 -10 0 10 20

The fourier transform of the hard-wall boundary condtion on the centre body is simply β/r(h,u)=0满足. To satisfy this boundary condition, we choose to write the solution of Bessel eqution as follows:

β(r,u)={A(u))Hm(1)(λ0wr) if r>1B(u)[Ym(λ1wh)Jm(λ1wr)Jm(λ1wh)Ym(λ1wr)] if h<r<1

where Jm and Ym are the Bessel and Neumann functions of order m, and Hm(1)=Jm+iYm is the Hankel function of the first type.

那么问题转化为,在复平面范围内,λ避开branch cut的情况下,如何保证λ0wr虚数部分始终为正数? 其中,w=wr+iwi=|w|eiε,0επ/2,即wr,wi均大于0

In [7]:
θ=range(0.0,stop=π/2-0.01,length=100)
plot(θ,tan.(θ))
Out[7]:
0.0 0.5 1.0 1.5 0 25 50 75 100 y1
In [8]:
wr = 0:1:100
wi = 0:1:100
f = z -> z
r = z -> abs(f(z))
θ = z -> angle(f(z))
plot( contourf(wr, wi, r.(wr' .+ im.*wi); title="abs"),
      contourf(wr, wi, θ.(wr' .+ im.*wi); title="phase"))
Out[8]:
0 20 40 60 80 100 0 20 40 60 80 100 abs 0 25 50 75 100 125 0 20 40 60 80 100 0 20 40 60 80 100 phase 0 0.25 0.50 0.75 1.00 1.25 1.50
In [52]:
M0=0 #M0=0,u±=±1
wr = 4
wi = 0
w=wr+ im*wi
phaseplot(-3..3, -3..3, z ->w*sqrt((1-z*M0)^2-z^2))
Out[52]:
-3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3
In [9]:
#u is in complex plane,for lanta_0,it has branch cuts
#λ0=sqrt((1-z*M0)^2-z^2)
#λ1=sqrt((C1-z*M0)^2-z^2)
M0=0 #M0=0,u±=±1
M1=0.5
C1=0.8
wr = 4
wi = 0
w=wr+ im*wi

xx = -10:0.01:10
yy = -10:0.01:10
f = z -> w*sqrt((1-z*M0)^2-z^2)

r = z -> abs(f(z))
θ = z -> angle(f(z))
image=z -> imag(f(z))
plot( contourf(xx, yy, r.(xx' .+ im.*yy); title="abs"),
      contourf(xx, yy, θ.(xx' .+ im.*yy); title="phase"),
      contourf(xx, yy, image.(xx' .+ im.*yy); title="image"))
Out[9]:
-10 -5 0 5 10 -10 -5 0 5 10 abs 0 10 20 30 40 50 -10 -5 0 5 10 -10 -5 0 5 10 phase - 1.5 - 1.0 - 0.5 0 0.5 1.0 1.5 -10 -5 0 5 10 -10 -5 0 5 10 image - 30 - 20 - 10 0 10 20 30
In [30]:
#u is in complex plane,for lanta_0,it has branch cuts
#λ0=sqrt((1-z*M0)^2-z^2)
#λ1=sqrt((C1-z*M0)^2-z^2)
M0=0.3 #M0=0,u±=±1
M1=0.5
C1=0.8
wr = 4
wi = 4
w=wr+ im*wi

xx = -10:0.01:10
yy = -10:0.01:10
f = z -> w*sqrt((1-z*M0)^2-z^2)
#r = z -> abs(f(z))
#θ = z -> angle(f(z))
image=z -> imag(f(z)) #求f函数的虚数部分

tanε=wi/wr
u1ip=range(10,stop=-10,length=10)
u1rp=u1ip./-tanε+ones(10,1)*1/(M0+1)
u1im=range(10,stop=-10,length=10)
u1rm=-u1im./tanε+ones(10,1)*1/(M0-1)
u2ip=range(10,stop=-10,length=10)
u2rp=u2ip./-tanε+ones(10,1)*C1/(C1*M0+1)
u2im=range(10,stop=-10,length=10)
u2rm=-u2im./tanε+ones(10,1)*1/(C1*M0-1)
plot(contourf(xx, yy, image.(xx' .+ im.*yy); title="image"))
plot!(u1rp,u1ip,aspect_ratio=:equal;label="R1+0",arrow=1)
plot!(u1rm,u1im,aspect_ratio=:equal;label="R1-0",arrow=1)

#plot!(u2rp,u2ip,aspect_ratio=:equal;label="R2+0")
#plot!(u2rm,u2im,aspect_ratio=:equal;label="R2-0")


#phaseplot(-3..3, -3..3, z ->w*sqrt((1-z*M0)^2-z^2))
Out[30]:
-10 -5 0 5 10 -10 -5 0 5 10 image - 25 0 25 50 75 R1+0 R1-0

Finally,we can derive the Winer-Hopf equation:

G(u)=G(u)+G+(u)=wK(u)F+(u)λ0λ1
where K is the Wiener-Hopf kernel:
K(u)=D1(1uM1)2λ0Ym(λ1wh)Jm(λ1w)Jm(λ1wh)Ym(λ1w)Ym(λ1wh)Jm(λ1w)Jm(λ1wh)Ym(λ1w)(1uM0)2λ1Hm(λ0w)Hm(1)(λ0w)
接下来,对K(u)的zeros和poles进行分析。

In [53]:
M0=0.25
M1=0.5
D1=C1=1
m=2
w=4+4im
h=0.1
λ0=(u)->sqrt((1-u*M0)^2-u^2)
λ1=(u)->sqrt((C1-u*M1)^2-u^2)
tanε=imag(w)/real(w)
u1ip=range(10,stop=-10,length=10)
u1rp=u1ip./-tanε+ones(10,1)*1/(M0+1)
u1im=range(10,stop=-10,length=10)
u1rm=-u1im./tanε+ones(10,1)*1/(M0-1)
u2ip=range(10,stop=-10,length=10)
u2rp=u2ip./-tanε+ones(10,1)*C1/(C1*M1+1)
u2im=range(10,stop=-10,length=10)
u2rm=-u2im./tanε+ones(10,1)*1/(C1*M1-1)

#####
𝐽𝑚_λ1w=(u) -> besselj(m, λ1(u)*w)
𝑌𝑚_λ1w=(u) -> bessely(m, λ1(u)*w)
#𝐽𝑚_λ0w=(u) -> besselj(m, λ0(u)*w)
#𝑌𝑚_λ0w=(u) -> bessely(m, λ0(u)*w)
𝐽′𝑚_λ1w=(u) -> 1/2*(besselj(m-1, λ1(u)*w)-besselj(m+1, λ1(u)*w))
𝑌′𝑚_λ1w=(u) -> 1/2*(bessely(m-1, λ1(u)*w)-bessely(m+1, λ1(u)*w))
𝐽′𝑚_λ1wh=(u) -> 1/2*(besselj(m-1, λ1(u)*w*h)-besselj(m+1, λ1(u)*w*h))
𝑌′𝑚_λ1wh=(u) -> 1/2*(bessely(m-1, λ1(u)*w*h)-bessely(m+1, λ1(u)*w*h))
𝐻1𝑚1_λ0w=(u) -> hankelh1(m, λ0(u)*w)
𝐻′1𝑚1_λ0w=(u) ->hankelh1(m-1, λ0(u)*w)-m/(λ0(u)*w)*hankelh1(m, λ0(u)*w)
#𝐻′1𝑚2_λ0w=(u) ->-hankelh1(m+1, λ0(u)*w)+m/(λ0(u)*w)*hankelh1(m, λ0(u)*w)
K=(u)->D1*(1-u*M1)^2*λ0(u)*(𝑌′𝑚_λ1wh(u)*𝐽𝑚_λ1w(u)-𝐽′𝑚_λ1wh(u)*𝑌𝑚_λ1w(u))/(𝑌′𝑚_λ1wh(u)*𝐽′𝑚_λ1w(u)-𝐽′𝑚_λ1wh(u)*𝑌′𝑚_λ1w(u))-(1-u*M0)^2*λ1(u)*𝐻1𝑚1_λ0w(u)/𝐻′1𝑚1_λ0w(u)
phaseplot(-15..15,-10..10, K)
plot!(u1rp,u1ip,aspect_ratio=:equal;label="R1+0",arrow=1,legend=:bottomright)
plot!(u1rm,u1im,aspect_ratio=:equal;label="R1-0",arrow=1)
plot!(u2rp,u2ip,aspect_ratio=:equal;label="R2+0",arrow=1)
plot!(u2rm,u2im,aspect_ratio=:equal;label="R2-0",arrow=1)
Out[53]:
-15 -10 -5 0 5 10 15 -10 -5 0 5 10 R1+0 R1-0 R2+0 R2-0